module points
   use precision
   private
   public operator(.dot.), sqrt, norm
   public operator(+), operator(-), operator(*), operator(/)

   type,public:: point
      real(adequate):: x = 0.0
	  real(adequate):: y = 0.0
	  real(adequate):: z = 0.0
   end type point

   interface operator(+)
      module procedure point_add, point_add_s, point_s_add
   end interface

   interface operator(-)
      module procedure point_subtract, point_subtract_s, &
	     point_s_subtract, point_negate
   end interface

   interface operator(*)
      module procedure point_multiply, point_multiply_s, &
	     point_s_multiply
   end interface

   interface operator(/)
      module procedure point_divide, point_divide_s, &
	     point_s_divide
   end interface

   interface norm
      module procedure point_norm
   end interface

   interface operator(.dot.)
      module procedure point_dot
   end interface

   interface sqrt
      module procedure point_sqrt
   end interface

contains

   function point_add (a, b) result(r)
   type(point), intent(in):: a, b
   type(point):: r

   r%x = a%x + b%x
   r%y = a%y + b%y
   r%z = a%z + b%z
   end function point_add

   function point_add_s (a, b) result(r)
   type(point), intent(in):: a
   real(adequate), intent(in):: b
   type(point):: r

   r%x = a%x + b
   r%y = a%y + b
   r%z = a%z + b
   end function point_add_s 

   function point_s_add (a, b) result(r)
   type(point), intent(in):: b
   real(adequate), intent(in):: a
   type(point):: r

   r%x = a + b%x
   r%y = a + b%y
   r%z = a + b%z
   end function point_s_add

   function point_negate (a) result(r)
   type(point), intent(in):: a
   type(point):: r

   r%x = -a%x
   r%y = -a%y
   r%z = -a%z
   end function point_negate

   function point_subtract (a, b) result(r)
   type(point), intent(in):: a, b
   type(point):: r

   r%x = a%x - b%x
   r%y = a%y - b%y
   r%z = a%z - b%z
   end function point_subtract

   function point_subtract_s (a, b) result(r)
   type(point), intent(in):: a
   real(adequate), intent(in):: b
   type(point):: r

   r%x = a%x - b
   r%y = a%y - b
   r%z = a%z - b
   end function point_subtract_s

   function point_s_subtract (a, b) result(r)
   type(point), intent(in):: b
   real(adequate), intent(in):: a
   type(point):: r

   r%x = a - b%x
   r%y = a - b%y
   r%z = a - b%z
   end function point_s_subtract


   function point_multiply (a, b) result(r)
   type(point), intent(in):: a, b
   type(point):: r

   r%x = a%x * b%x
   r%y = a%y * b%y
   r%z = a%z * b%z
   end function point_multiply

   function point_multiply_s (a, b) result(r)
   type(point), intent(in):: a
   real(adequate), intent(in):: b
   type(point):: r

   r%x = a%x * b
   r%y = a%y * b
   r%z = a%z * b
   end function point_multiply_s

   function point_s_multiply (a, b) result(r)
   type(point), intent(in):: b
   real(adequate), intent(in):: a
   type(point):: r

   r%x = a * b%x
   r%y = a * b%y
   r%z = a * b%z
   end function point_s_multiply

   function point_divide (a, b) result(r)
   type(point), intent(in):: a, b
   type(point):: r

   r%x = a%x / b%x
   r%y = a%y / b%y
   r%z = a%z / b%z
   end function point_divide

   function point_divide_s (a, b) result(r)
   type(point), intent(in):: a
   real(adequate), intent(in):: b
   type(point):: r

   r%x = a%x / b
   r%y = a%y / b
   r%z = a%z / b
   end function point_divide_s

   function point_s_divide (a, b) result(r)
   type(point), intent(in):: b
   real(adequate), intent(in):: a
   type(point):: r

   r%x = a / b%x
   r%y = a / b%y
   r%z = a / b%z
   end function point_s_divide

   function point_norm(self) result(r)
   type(point),intent(in):: self
   real(adequate) r
   r = sqrt(self%x**2 + self%y**2 + self%z**2)
   end function point_norm

   function point_dot(a, b) result(r)
   type(point), intent(in):: a, b
   real(adequate):: r
   real(precise):: s
   s = real(a%x, precise) * b%x + real(a%y,precise) * b%y + &
       real(a%z, precise) * b%z
   r = real(s, adequate)
   end function point_dot

   function point_sqrt(a) result(r)
   type(point), intent(in):: a
   type(point):: r
   r%x = sqrt(a%x)
   r%y = sqrt(a%y)
   r%z = sqrt(a%z)
   end function point_sqrt
end module points